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Abstract 

We revisit the Mittag-Leffler functions of a real variable t, with one, 
two and three order-parameters {a,/?, 7}, as far as their Laplace 
transform pairs and complete monotonicty properties are concerned. 
These functions, subjected to the requirement to be completely 
monotone for t > 0, are shown to be suitable models for non- 
Debye relaxation phenomena in dielectrics including as particular 
cases the classical models referred to as Cole-Cole, Davidson-Cole 
and Havriliak-Negami. We show 3D plots of the response functions 
and of the corresponding spectral distributions, keeping fixed one of 
the three order-parameters. 
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Foreword 



This E-print reproduces the revised version of the paper published in EPJ- 
ST, Vol. 193 (2011), pp. 161-171. The revision concerns the proper use 
of the terms relaxation function and response function in the literature on 
dielectrics. In the published paper, starting from Eq (1.1), the authors had 
referred the inverse Laplace transform of the complex susceptibility as the 
relaxation function. This is not correct because the inversion provides the 
so-called response function as pointed out to FM by Prof. Karina WERON 
(KW) to whom the authors are very grateful. As a matter of fact the 
relationship between the response function and the relaxation function can 
better be clarified by their probabilistic interpretation investigated in several 
papers by KW. As a consequence, interpreting the relaxation function as a 
survival probability ^(t), the response function turns out to be the probability 
density function corresponding to the cumulative probability function <E>(t) = 
1 — ty(t). Then, denoting by 4>(t) the response function, we have 

As KW has pointed out, both functions have very different properties and 
describe different physical magnitudes; only in the Debye (pure exponential) 
case the properties coincide. The relaxation function describes the decay of 
polarization whereas the response function its decay rate (the depolarization 
current). However, for physical realizability, both functions are required to 
be completely monotone with a proper spectral distribution so our analysis 
can be properly transferred from response functions to the corresponding 
relaxation functions, whereas the corresponding cumulative probability 
functions turn out to be Bernstein (or creep) functions, that is positive 
functions with a completely monotone derivative. 



1 Introduction 

It is well recognized that relaxation phenomena in dielectrics deviate more 
or less strongly from the classical Debye law for which the Laplace transform 
pair for complex susceptibility (s = —iui) and response function it > 0) reads 
in an obvious notation, 

1^ -5-e D (*) = e _ *. (1.1) 
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Here, for the sake of simplicity, we have assumed the frequency co and 
the time t normalized with respect to a characteristic frequency oud and 
a corresponding relaxation time Td = 1/ood- 

In the literature a number of laws have been proposed to describe the non- 
Debye (or anomalous) relaxation phenomena in dielectrics, of which the most 
relevant ones are referred to Cole - Cole (C-C), Davidson - Cole (D-C) and 
Havriliak - Negami (H-N) laws , see e.g., the classical books by Jonscher 
[13, 14]. Several authors have investigated these laws from different points 
of view, including Karina Weron and her associates, see e.g., [15, 16, 29, 30], 
Hilfer [11, 12] and Hanyga and Seredyhska [7]. 

In particular, Hilfer surveyed the analytical expressions in the frequency and 
time domain for the main non-Debye relaxation processes and provided the 
response functions corresponding to the complex frequency-dependent Cole- 
Cole, Davidson-Cole and Havriliak-Negami susceptibilities in terms of Fox 
if— functions. This class of functions is quite general so it includes the 
Mittag-Leffler functions that we prefer to use to characterize the above laws 
in a more accessible way. 

On the other hand, for linear systems, the connection between weak 
dissipativity and positive defmiteness of the response functions as well as 
between monotone energy decay and complete monotonicity of the response 
functions were discussed by Hanyga and Seredyhska [7], and in references 
therein, in terms of functions of Mittag-Leffler type. A subordination model 
of anomalous diffusion leading to the two-power-law relaxation responses has 
been proposed by Stanislavski et al. [30], where the authors have presented 
a novel two-power relaxation law and shown its relationship to the H-N law 
by using functions of Mittag-Leffler type. 

Moreover, Novikov et al. [23] have discussed anomalous relaxation in 
dielectrics providing a differential equation with fractional derivatives to 
describe the relaxations of the C-C, C-D and H-N types in dielectrics. The 
corresponding solutions were presented in terms of Mittag-Leffler functions. 

Fractional relaxation and governing equations for dielectrics characterized 
by the H-N response function has been pointed out by Sibatov and Uchaikin 
[28] who have presented a comparison of numerical results with experimental 
data. 

The purpose of this paper is to present a general model for anomalous 
relaxation in dielectrics that includes as particular cases the classical C-C, 
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D-C and N-H laws. Our model is still based on the Mittag-Lefner functions 
but depending on three order-parameters that, in view of their complete 
monotonicity in the time domain, ensure the existence of a suitable spectrum 
of relaxation times required for the physical realizability. On the realizability 
requirements the reader can be addressed, in addition to [7], to the treatise 
by Zemanian [31]. 

In a subsequent paper we will investigate how to provide the differential 
equation of fractional order that is expected to govern our general model, so 
generalizing the approach and the results by Novikov et al. 

The present paper is organized as follows. 

In section 2 we recall the definitions of the most common functions of Mittag- 
Leffler type, namely those depending on one, two and three order-parameters 
{a, (3, 7}. For these functions we exhibit the corresponding Laplace transform 
pairs and we discuss their properties of complete monotonicity. We show 
that the complete monotonicity is ensured if the independent variable is real 
and negative and the three order-parameters are subjected to the conditions 
< a,/3,7 < 1 with 0:7 < (5. For other properties of the Mittag-Leflier 
functions we refer the reader mainly to texts on Special Functions and 
Fractional Calculus, including i.e. [4, 10, 17, 18, 20, 24], just to cite some 
of the most recent ones. In fact, as pointed out by Gorenflo and Mainardi 
[5, 19], functions of Mittag-Leffler type enter as solutions of many problems 
dealt with fractional calculus so that they like to refer to the Mittag-Leflier 
function to as the Queen function of Fractional Calculus, in contrast with its 
role of a Cinderella function played in the past. 

In section 3 we show that for special cases of the triplet {a, (3, 7} with 
a 7 = (3 our functions provide the response functions for the classical models 
of Cole-Cole, Davidson-Cole, and Havriliak-Negami. As a consequence, we 
expect that in the more general case 0:7 < f3 the corresponding Mittag- 
Leffler functions, being completely monotone, can provide further models 
for processes of anomalous (non-Debye) relaxation in dielectrics. For some 
study-cases, taking fixed two of the three order-parameters, we provide 
3D plots of the responses functions and of the corresponding spectral 
distributions, in order to better visualize the positivity and the variability of 
the considered functions. 

Finally, section 4 is devoted to conclusions and final remarks together with 
a discussion on the direction of future work. 
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2 The Mittag-Leffler functions 



2.1 Definitions 

We start to recall the definition in the complex plane of the generalized 
Mittag-Leffler function introduced by Prabhakar [26], known as Prabhakar 
function or 3-parameter Mittag-Leffler function 

E IA Z ) ■= E ^, r /Zh m zn . Re M > 0» Re W > 0. Re M > 0, (2.1) 

n=0 ^ 

where 

( 7 ) n = 7 ( 7 + l)...( 7 + n-l) = ^±^. 

It turns out to be an entire function. For 7 = 1 we recover the 2-parameter 
Mittag-Leffler function (also known as Wiman or Agarwal generalized 
Mittag-Leffler function) 

00 n 

(2 ' 2) 

and for 7 = j3 = 1 we recover the standard Mittag-Leffler function 

EJz) :=Y -. (2.3) 

v ; ^f«n+l 



n=0 



2.2 The Laplace transform pairs 

Let us now consider the relevant formulas of Laplace transform pairs related 
to the above three functions, already known in the literature for a, (5 > 
and < 7 < 1, when the independent variable is real of type at where t > 
is interpreted as time and Qj dbS db certain constant. 

Let us start with the most general function. Substituting the series 
representation of the Prabhakar generalized Mittag-Leffler function in the 
Laplace transformation yields the identity 

e-^-^nd^E^G)''. (2.4) 

71=0 ^ 1 ' 
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On the other hand (binomial series) 

(i + *r = f w r(1 " 7 \ , * B = Yi- 1 ) 11 ^STr - n • (2-5) 

n=0 V ' y n=0 V ' ; 

Comparison of Eq. (2.4) and Eq. (2.5) yields the Laplace transform pair 

tP - 1 El 3 (at a ) + — ^ = ^ . (2.6) 

Eq. (2.6) holds (by analytic continuation) for Re(a) > 0, Re(/3) > 0. 

In particular we get the known Laplace transform pairs, see e.g. [18, 24], 

tP ~ 1 £ a fl(at a ) -T- = — ^ , (2.7) 

' py ' s a -a 1 - as~ a ' v 7 

„a— 1 „— 1 

£ a (at Q ) - = . (2.8) 



2.3 Complete monotonicity 

We recall that a function fit) with t > is completely monotone (CM) if it 
is positive and its derivatives are alternating in sign, namely 

(-l)"/ (n) (0 > 0, t > 0. (2.9) 

Thus, for the Bernstein theorem that states a necessary and sufficient 
condition for the CM, the function can be expressed as a real Laplace 
transform of non-negative (generalized) function, namely 

poo 

f(t)= e~ rt K(r) dr , K(r) > , t > . (2.10) 
Jo 

By the way, the determination of such non-negative function K{r) (the 
Laplace measure) is a standard method to prove the CM of a given function 
defined in the positive real axis 1R + . In physical applications the function 
K(r) is usually referred to as the spectral distribution function, in that it 
is related to the fact that the process governed by f(t) can be expressed 
in terms of a continuous distribution of elementary (exponential) relaxation 
processes. 
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For the Mittag-Leffler functions in one and two-order parameters the 
conditions to be CM were proved respectively by Pollard [25] in 1948 and by 
Schneider [27] in 1996, see also Miller and Samko [21, 22] for further details. 
These conditions require the independent variable to be real and negative 
and the order-parameters such that < a < 1 and a < (3 < 1, respectively. 

A method for proving CM suitable for our Mittag-Leffler functions (2.6) 
(2.7) (2.8) is that to arrive at the non-negative spectral distribution K(r) 
by inverting the corresponding Laplace transform by using the Bromwich 
countour integral. Indeed this was the method followed by Gorenflo and 
Mainardi [5] to prove the CM of E a (—t~ a ) for < a < 1 and to determine 
the corresponding spectral function. We recall the result 

E (-t a ) = / V< K (r) dr K (r) = — > 

(2.11) 

As a matter of fact the function K a (r) is derived as an exercise in complex 
analysis by evaluating the contribution on the branch cut (the negative real 
axis) of the Bromwich integral and turns out to be provided by the so-called 
Titchmarsch formula 



1 f s Q_1 
K a (r) = —lm\ —- 

7T I S a + 1 



(2.12) 



We now prove that the function 

fo(0:=^ _1 ^(-O, (2-13) 
with Laplace transform (as derived from (2.7) with a = — 1) 

&W = ( S JTTF- < 2 ' 14 > 

where the notation £g(0 been introduced for future convenience, is CM iff 

0<a,/3,7<l with «7</3. (2.15) 

Furthermore, by using the Titchmarsh formula we will derive its spectral 
distribution function. For this purpose we take advantage of the requirements 
stated in the treatise by Gripenberg et al. [6], see Theorem 2.6, pp. 144- 
145, that provide necessary and sufficient conditions to ensure the CM of a 
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function f(t) based on its Laplace transform f(s). Hereafter, we recall this 
theorem by using our notation. 

Theorem The Laplace transform f(s) of a function f(t) that is locally 
integrable on 1R + and CM the following properties: 

(i) f(s) an analytical extension to the region C — IRT; 

(ii) J{x) = f*(x) forxe (0,oo); 

(iii) lim f(x) = 0; 

x— >oo 

(iv) Im{/(s)} < for Re{s} > 0; 

(v) Im{s f(s)} > for Im{s} > and f(x) > for x G (0, oo). 
Conversely, every function f(s) that satisfies (i)-(iii) together with (iv) or 
(v), is the Laplace transform of a function f(t), which is locally integrable 
on 1R + and CM on (0,oo). 

To prove the CM of £g(0 we must recognize that all the requirements (i)-(v) 
are satisfied for £g( s ) under the conditions (2.15). We first recognize that for 

< a < 1, < p < 1, < 7 < 1, 

the requirements (i)-(ii) are satisfied. The requirement (iii) is satisfied for 
the additional condition — (3 < 0. Then it remains to prove that the 
requirement (v) is satisfied for the function 

s£ G (s) = - =s 1 -P{- . (2.16) 

V ' (S a + 1)7 \1 + S~ a J K ' 

Putting s = re'* with < < tc, upper half C + (Im{s} > 0) we can get 

s(o(s) = ^-k Q+I-LJ . (2.17) 
Taking its imaginary part we can write 

Im{s^ G (s)} = ^Im{(re^) 1 - /3 (1 + r~ a e i(f)a y} , with Q = |l + s~ a | 27 >0. 

(2-18) 

Thus, we must proof that this imaginary part is always non- negative in the 
upper half-plane C + . Setting 

Zl : = r 1 "^ 1 "^ , z 2 :=r-V^, 



8 



we recognize that Z\ lies in C + being < (5 < 1, that z 2 and consequently 
1 + ^2 lies in C + being < a < 1 and the same for (1 + z-if 1 being < 7 < 1. 
We therefore conclude that under these requirement on {a (3, 7} we get 

Im{s£ G (s)} > 0. (2.19) 

As a matter of fact we have proved a noteworthy result stated without explicit 
prove by Mainardi [18] in his recent book. For this purpose we have adopted 
the same method followed by Hanyga and Seredyhska [7] who have limited 
themselves to prove, starting form the Laplace transform pair (2.6) with 
f3 — 1, and basing on [6], 

E l,i(~t a ) CM if 0<«<1,0<7<1. (2.20) 
So we have extended their result for < j3 < 1. 

2.4 The spectral distribution function 

The Bernstein theorem ensures that our response function £g(^)> being CM, 
can be expressed with a spectral distribution function K^^r) > such that 

POO 

&(*)= / e- rt K^(r)dr. (2.21) 
Jo 

In view of Titchmarsh formula (2.12) applied to the Laplace transform (2.14) 
we get 

Kl Jr) := — Im <^ } . (2.22) 

We note that for (3 = 7 = 1 we obtain from (2.22) the spectral distribution 
function of the classical Mittag-Leffler function outlined in (2.11). 
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3 Mathematical models for dielectric relax- 
ation 



We intend to discuss how the Laplace transform pair outlined in Eqs. (2.13)- 
(2.14) coupled with the condition (2.15) on the 3-order parameters can be 
considered as the pair £g(£) (t > 0) and £g( s ) ( s — iw) for a possible 
mathematical model of the response function and the complex susceptibility 
in the framework of a general relaxation theory of dielectrics. We first show 
how the three classical models referred to Cole-Cole (C-C), Davidson-Cole 
(D-C) and Havriliak-Negami (H-N) for cry = /3 are contained in our general 
model when cry = (3 according to the scheme 

[ 0<a<l )J 9 = a,7 = l C-C {a} , 
cry = /3 with < a = 1 , /3 = 7 , < 7 < 1 D-C {7} , (3.1) 
[o<a<l,0<7<l H-N {a, 7} . 

Then, we consider some study-cases when the inequality cry < (3 holds. For 
this purpose we exhibit 3D plots for the response function £g(0> see (2-13), 
and the corresponding spectral distribution K^^r), see (2.21) keeping fixed 
two of the three-order parameters. 

3.1 The classical dielectric functions 

The Cole— Cole relaxation model The C-C relaxation model is a non- 
Debye relaxation model depending on one parameter, say a (0 < a < 1), see 
[1, 2], that for a — 1 reduces to the standard Debye model. 

ec-coo = - £c-c(0 = t*- 1 Kj-n = -^-(-o , o < a < 1 . 

(3.2) 

The Davidson— Cole relaxation model The D-C relaxation model is a 
non-Debye relaxation model depending on one parameter, say 7 (0 < 7 < 1), 
see [3], that for 7 = 1 reduces to the standard Debye model. The 
corresponding complex susceptibility (s = —iu) and response function read 

&-C(*) = (T^yr ^ *D-C(0 = ^ E IH) = e " ' ° < 7 < 1 • 

(3.3) 
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The Havriliak— Negami relaxation model The H-N relaxation model 
is a non-Debye relaxation model depending on two parameters, say a (0 < 
a < 1) and 7 (0 < 7 < 1), see [8, 9], that for a = 7 = 1 reduces to 
the standard Debye model. The corresponding complex suceptibility and 
response function read 

*H-N(*) = (TT^f + *H-N(0 = ^ ^, Q7 (-0 , < a, 7 < 1 . (3.4) 

We recognize that the H-N relaxation model for 7 = 1 reduces to the C- 
C model, while for a — 1 to the D-C model. We also note that whereas 
for the C-C and H-N models the corresponding response functions decay 
like a certain negative power of time, the D-C response function exhibits an 
exponential decay. 

3.2 Survey of the general response functions with their 
spectral distributions 

In order to visualize the effects of varying the three order-parameters in 
our general model we survey some particular cases by exhibiting separately 
the response functions £g(*) for < t < 10 given by Eq. (2.21) and the 
corresponding spectral distributions for < r < 10 given by Eq. (2.22). 
For this purpose we provide 3D plots by keeping fixed two of the three 
order-parameters {a,/?, 7}, all of them less than unity and subjected to the 
condition a 7 < (5. 

In Figs. 1. and 2., for fixed a = 1/2 we compare versus 7 G (0, 1) the plots 
of £g(*) an d of the corresponding spectral distribution for /3 — 7/2, 27/3. 
In Figs. 3. and 4., for fixed 7 = 1/2 we compare versus a G (0, 1) the plots 
of £g(£) an d of the corresponding spectral distribution for f3 = a/2, 2a/3. 

In Figs. 5. and 6., for fixed a = 1/2,7 — 1 an d a = 2/3,7 = 1/2 we 
compare versus (3 G (1/3,1) the plots of £g(£) an d of the corresponding 
spectral distribution. 
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Figure 1: Plot of the response Function £g(£) versus 7 G (0, 1) and t G (0, 10) 
for {a = 1/2,(3 = 7/2,7} (dark grey/blue surface) and {a = 1/2,(3 = 
27/3,7} (light grey/yellow surface). 
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Figure 2: Plot of the Spectral Distribution Function Ky 2 ^^ 2 (r) (dark 
grey/blue surface) and K^ 22 ^^ 3 (r) (light grey/yellow surface) for < 7 < 1 
and < r < 10. 
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0.0 




Figure 3: Plot of the response Function £c?(i) versus a G (0, 1) and t G (0, 10) 
for {a, (3 = a/2, 7 = 1/2} (dark grey/blue surface) and {a, f3 = 2a/3, 7 = 
1/2} (light grey /yellow surface). 
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Figure 4: Plot of the Spectral Distribution Function Kj a j 2 {r) (dark 

1 /2 

grey/blue surface) and K a ' 2a ^ 3 (r) (light grey/yellow surface) for < a < 1 
and < r < 10. 
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Figure 5: Plot of the response Function £g>(£) versus G (1/3,1) and t G 
(0, 10) for {a = 1/2, (3, 7 = 1} (dark grey/blue surface) and {a = 2/3,(3 = 
2a/3, 7 = 1/2} (light grey /yellow surface). 
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Figure 6: Plot of the Spectral Distribution Function K\j 2 ^(r) (dark 

1/2 

grey/blue surface) and K 2 L^{r) (light grey/yellow surface) for 1/3 < (3 < 1 

1 /2 

and < r < 10. Note that K 2 ', 3 ^(r) < when the condition /3 > 07 is not 
satisfied. 
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4 Conclusions 



In this paper we have presented a quite general mathematical model for 
dielectrics that exhibit deviations from the standard relaxation Debye law. 
This model is based on a time response function expressed in terms of a 
Mittag-Leffler function with three order-parameters. A restriction on these 
parameters is required to ensure its complete monotonicity for t > 0, so that 
the resulting relaxation process can be seen as a continuous distribution of 
elementary exponential processes by means of a corresponding (non-negative) 
spectral distribution function. 

Our approach allows one to derive from a unique mathematical framework 
the classical Cole-Cole, Davidson-Cole and Havriliak-Negami models that 
are usually adopted in the literature. However, other laws can be derived 
from our model that could better fit some experimental data. 

For some study-cases we have exhibited plots of the response functions 
(versus t) along with their corresponding spectral distributions (versus r) 
keeping fixed two of the three order-parameters {a, (3, 7} of our Mittag- 
Leffler function. We can note that when the condition < «7 < /3 < 1 
is not satisfied the spectral distribution functions turn out to be negative. 

In a future paper we are interested to provide the fractional differential 
equation governing our general response function extending the results 
obtained by Novikov et al. [23] and by Stanislavski et al. [30] for the 
Havriliak-Negami model. 
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